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ABSTRACT.  In  October  2010  the  US  Naval  Observatory  together  with  the  Space 
Research  Centre  in  Warsaw  initiated  the  Earth  Orientation  Parameters  Combination  of 
Prediction  Pilot  Project,  which  was  accepted  by  the  IERS  Directing  Board.  The  goal  of  this 
project  is  to  determine  the  feasibility  of  combining  Earth  Orientation  Parameters  (EOP) 
predictions  on  an  operational  basis.  The  ensemble  predictions  of  EOPs  are  more  accurate 
than  the  results  from  individual  predictions.  The  pole  coordinate  data  predictions  from 
different  prediction  contributors  and  ensemble  predictions  computed  by  the  U.S.  Naval 
Observatory  are  studied  to  determine  the  statistical  properties  of  polar  motion  forecasts  by 
looking  at  second,  third  and  fourth  moments  about  the  mean.  The  increase  of  prediction  errors 
in  pole  coordinate  data  can  be  due  to  the  change  of  phase  of  the  annual  oscillation  in  the  joint 
atmospheric-ocean  excitation  function. 
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1.  INTRODUCTION 

The  International  Celestial  and  Terrestrial  Reference  Systems  (ITRS  and  ICRS)  are 
defined  by  coordinates  of  extragalactic  radio  sources  observed  by  VLBI  and  by  the 
coordinates  and  velocities  of  sites,  respectively,  determined  from  space  geodetic  techniques 
e.g.  GNSS,  DORIS  (Doppler  Orbitography  and  Radiopositioning  Integrated  by  Satellite), 
SLR,  LLR  and  VLBI.  The  EOP  predictions  enable  the  real  time  transformation  between  the 
realizations  of  these  systems  that  are  the  International  Celestial  and  Terrestrial  Reference 
Frames  (ICRF  and  ITRF).  During  the  last  decades  the  definition  and  practical  realization  of 
these  systems  has  been  adapting  to  increasing  accuracy  of  the  measurements  provided  by 
these  techniques.  The  real  time  transformation  between  these  two  frames  is  a  function  of  time 
realized  by  predictions  of  x,  y  pole  coordinates,  UT1-UTC  and  a  precesion-nutation 
extrapolation  model  together  with  the  predictions  of  dX,  dY  precesion-nutation  residuals.  The 
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Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 
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determination  accuracy  of  these  EOPs  observed  by  space  geodetic  techniques  is  of  the  order 
of  few  millimeters  on  the  Earth’s  surface,  however  such  accuracy  is  not  possible  in  the  real 
time  transformation  due  to  the  latency  of  products  provided  by  these  techniques.  Thus,  it  is 
necessary  to  predict  these  transformation  parameters.  However  their  prediction  error  even  for 
a  few  days  in  the  future  is  several  times  less  than  their  determination  accuracy  (Kosek  at  al. 
2004,  Luzum  2010).  Predictions  of  EOPs  are  now  provided  by  the  IERS  Rapid  Service 
Prediction  Centre  in  USNO,  Washington  DC.  The  x,  y  pole  coordinates  are  predicted  using  a 
combination  of  LS  extrapolation  and  autoregressive  prediction  of  the  least-squares  (LS) 
residuals  (Kosek  et  al.  2004).  The  short-term  UTI  -UTC  data  are  predicted  using  the  forecast 
of  the  axial  component  of  atmospheric  angular  momentum  obtained  during  the  estimation  of  a 
global  dynamic  atmospheric  circulation  model  (Johnson  et  al.  2005).  The  accuracy  of  the  new 
IAU2000  precession-nutation  model  is  very  high,  thus  its  differences  with  respect  to 
observations  are  very  small  and  show  a  quasi-periodic  signal  related  to  the  free  core  nutation 
(FCN)  with  an  amplitude  not  exceeding  0.3  mas.  The  accuracy  of  the  real  time  transformation 
between  the  celestial  and  terrestrial  reference  frames  mostly  depends  on  the  prediction 
accuracy  of  the  x,  y  pole  coordinates  and  UT1-UTC  data. 

The  predictions  of  EOPs  are  being  used  by  the  NASA  Deep  Space  Network  (DSN) 
which  consists  of  the  network  of  antennas  located  in  1)  Goldstone,  California,  Mojave  desert; 
2)  near  Madrid,  Spain;  and  3)  near  Canberra,  Australia  to  control  interplanetary  missions, 
radio  and  radar  astronomical  observations  as  well  as  some  Earth  orbiting  satellites.  The 
longitude  distance  of  about  120°  between  these  sites  enables  observations  of  space  objects 
regardless  of  Earth  rotation.  High  accuracy  prediction  of  EOPs  becomes  very  important  for 
the  realization  of  altimetric  and  gravimetric  missions  as  well  as  the  determination  of  precise 
directions  of  extragalactic  radio  sources  observed  by  VLBI. 

2.  INTERNATIONAL  COOPERATION 

To  improve  the  EOP  prediction  accuracy,  international  cooperation  has  become 
necessary  due  to  degradation  of  the  EOP  prediction  accuracy  with  prediction  length.  This 
prediction  accuracy  degradation  is  mostly  caused  by  variable  short  period  oscillations  in  x,  y 
pole  coordinates  and  UTI  -UTC  data.  Such  wideband,  short-period  oscillations  in  pole 
coordinates  and  universal  time  data  are  excited  by  wideband  high-frequency  oscillations  of 
equatorial  components  of  the  joint  atmospheric-ocean  excitation  function  and  the  axial 
component  of  atmospheric  excitation  function,  respectively. 

The  first  campaign,  named  Earth  Orientation  Parameters  Prediction  Comparison 
Campaign  (EOPPCC),  was  initiated  at  the  Technical  University  of  Vienna  by  Prof.  H.  Schuh 
and  Space  Research  Centre  of  Polish  Academy  of  Sciences  (SRC  PAS)  by  W.  Kosek  and  M. 
Kalarus.  The  goal  of  this  campaign,  which  started  in  Oct.  2005  and  terminated  in  Mar.  2008, 
was  the  comparison  of  the  EOP  prediction  results  computed  by  different  prediction 
techniques  using  different  input  data.  About  10  participants  of  this  campaign  computed  EOP 
predictions  every  week  that  were  then  analyzed  at  the  SRC  PAS.  The  most  important 
achievement  of  this  campaign  was  to  show  that  ensemble  predictions  of  x,  y  pole  coordinate 
data  are  more  accurate  than  the  individual  ones  and  the  most  accurate  short-term  prediction 
method  for  UTI -UTC  was  the  Kalman  filter  which  used  the  10-day  forecast  of  the  axial 
component  of  atmospheric  angular  momentum  as  an  input  together  with  UTI -UTC  data 
(Kalarus  et  al.  2010). 

In  Apr.  2006  at  the  EGU  the  IERS  Working  Group  on  Predictions  (WGP)  was  initiated  to 
study  the  advantages  and  shortfalls  of  different  EOP  prediction  algorithms  as  well  as  to 
determine  what  prediction  products  are  useful  to  the  user  community  in  addition  to  making  a 
detailed  examination  of  the  fundamental  properties  of  the  different  input  data  sets.  Due  to  the 
initiative  of  the  WGP  members,  in  Oct.  2009  the  IERS  Workshop  on  Combination  and 
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Prediction  was  held  in  Warsaw.  The  goals  of  this  Workshop  were  to  determine  the  current 
state-of-the-art  for  the  EOP  prediction  in  terms  of  data  sets  and  algorithms  as  well  as  to 
discuss  recommended  actions  for  improving  EOP  predictions,  including  improving  the  IERS 
rapid  combined  series.  This  Workshop  generated  about  20  recommendations  related  to 
observations,  analysis  and  prediction  of  EOPs.  Following  these  recommendations  in  Oct. 
2010  the  Earth  Orientation  Parameters  Combination  of  Prediction  Pilot  Project  (EOPCPPP), 
was  initiated  by  the  IERS  Directing  Board.  The  goal  of  this  Project  is  to  determine  the 
feasibility  and  benefits  of  combining  EOP  predictions  on  a  daily  basis  and  to  determine  the 
best  algorithms  for  EOP  prediction  combinations  (USNO  2011,  SRC  2011).  In  this  Project 
about  10  participants  provide  EOP  predictions  ranging  from  roughly  90  days  up  to  one  year 
(Tab.  1). 


Table  1.  The  participants  of  the  EOPCPPP  and  their  contribution.  See 
http://maia.usno.navy.mil/eopcppp/  or  http://www.cbk.waw.pl/eopcppp  for  additional  details. 
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3.  ANALYSIS 

To  analyze  prediction  results  of  the  individual  participants  of  the  project  as  well  as 
ensemble  predictions  computed  daily  at  the  US  Naval  Observatory  the  following  statistics 
were  applied:  standard  deviation  (SD),  mean  absolute  error  (MEA),  skewness  (SKE)  and 
kurtosis  (KUR).  The  prediction  accuracy  can  be  expressed  using  the  first  two  statistics:  SD 
and  MAE.  The  SD  of  predictions  computed  by  the  individual  participant  or  for  the  ensemble 
prediction  for  the  ith  day  in  the  future  is  computed  by  the  following  formula: 

SD,  =  Jl /np  kT  - )2  >  i  =  (1) 

where  jul  j  =  1 jnp  {x[Jd  ~  x°bs ),  i  =  1,2,.. .M  ,  is  the  first  moment  about  the  mean,  x-"Jd 

is  the  jth  prediction  computed  for  the  i,h  day  in  the  future,  xobs  is  the  observed  future  value  of 
time  series,  M  is  the  maximum  prediction  length,  np  is  the  number  of  predictions. 

The  statistical  error  of  SD  for  the  ith  day  in  the  future  is  computed  by  the  following 
formula: 


&{SD),=SD,lfiTr. 


(2) 


The  MAE  of  predictions  for  the  ith  day  in  the  future  is  computed  by  the  following 
formula: 


MAE,  =  1 Jnp  X":,  <jd -x°bS  >  i  =  1,2,..., M . 


(3) 


The  statistical  error  of  MAE  for  the  ih  day  in  the  future  is  computed  by  the  following 
formula  (Kalarus  et  al.  2010): 


(j(MAE),  =  SD^{K-2)l(K-np). 


(4) 


Skewness  (SKE)  is  a  measure  of  the  asymmetry  of  the  probability  distribution  of  a  real¬ 
valued  random  variable.  Negative  skewness  indicates  that  the  tail  on  the  left  side  of  the 
probability  density  function  is  longer  than  the  right  side.  If  the  distribution  is  symmetric  then 
skewness  is  equal  to  zero.  The  SKE  of  predictions  for  the  ith  day  in  the  future  are  computed 
by  the  following  formula: 


SKE,  =  E 
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The  statistical  error  of  SKE  for  the  ih  day  in  the  future  depends  only  on  the  number  of 
predictions  and  is  computed  by  the  following  formula: 


p  ' 


&{SKE)i  =  yjsJn 


(6) 
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Kurtosis  (KUR)  or  excess  kurtosis  is  a  measure  of  the  "peakedness"  of  the  probability 
distribution  of  a  real-valued  random  variable.  If  kurtosis  is  equal  to  0  then  the  differences 
between  pole  coordinate  data  and  their  predictions  follow  a  normal  distribution.  If  it  is 
negative  then  the  probability  distribution  is  more  flat  which  means  that  there  are  a  larger 
number  of  extreme  deviations  in  the  variance  and  if  its  positive  then  the  number  of  large 
deviations  in  the  variance  is  less  than  in  the  case  of  normal  distribution.  The  kurtosis  of 
predictions  for  the  ith  day  in  the  future  are  computed  by  the  following  formula: 


KUR,  =  E 
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■3.  (7) 


The  statistical  error  of  KUR  for  the  ith  day  in  the  future  is  computed  by  the  following 
formula: 


&(KUR)i  =  y]24/  ~np  . 


(8) 


n  in  nn  on  itn  cn  zn  on  on  nn 


Day  s  in  the  future  Day  s  in  the  future 

Fig.  1.  Mean  absolute  error,  standard  deviation,  skewness  and  kurtosis  together  with  their 
error  bars  of  x  (blue),  y  (red)  predictions  computed  by  WK. 
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Fig.  2.  Mean  absolute  error,  standard  deviation,  skewness  and  kurtosis  together  with  their 
error  bars  of  x  (blue),  y  (red)  predictions  computed  by  BL. 
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Day  s  in  the  future  Day  s  in  the  future 

Fig.  3.  Mean  absolute  error,  standard  deviation,  skewness  and  kurtosis  together  with  their 
error  bars  of  x  (blue),  y  (red)  predictions  computed  by  VT. 
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Fig.  4.  Mean  absolute  error,  standard  deviation,  skewness  and  kurtosis  together  with  their 
error  bars  of  x  (blue),  y  (red)  predictions  computed  by  ZM. 


arcsec 


Days  in  the  future 


Days  in  the  future 


arcsec 


Days  in  the  future 


Days  in  the  future 


Fig.  5.  Mean  absolute  error,  standard  deviation,  skewness  and  kurtosis  together  with  their 
error  bars  of  x  (blue),  y  (red)  ensemble  predictions  computed  by  BL-ep. 
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All  four  statistics  were  computed  for  x,  y  pole  coordinate  data  predictions  provided  by 
BL,  WK,  VT  and  ZM  as  well  as  ensemble  predictions  of  these  data  provided  by  BL-ep. 
Figures  1-5  show  the  values  of  SD,  MAE,  SKE  and  KUR  as  a  function  of  prediction  length 
for  chosen  participants  of  the  project.  It  can  be  seen  that  for  different  participants/prediction 
methods  the  skewness  changes  are  similar.  The  negative  value  of  skewness  for  the  x  pole 
coordinate  is  associated  with  a  positive  value  for  the  y  pole  coordinate  except  for  participant 
ZM  where  the  skewness  values  are  negative  (Fig.  4).  Negative  or  positive  values  of  skewness 
suggest  that  the  probability  distribution  of  the  differences  between  the  pole  coordinate  data 
and  their  predictions  are  not  symmetrical.  A  similar  change  of  skewness  suggests  also  that 
there  are  irregular  variations  in  the  pole  coordinate  data  that  are  difficult  to  predict  by  any 
method. 

The  kurtosis  values  for  different  participants  are  usually  similar  and  decrease  with  the 
prediction  length  except  the  y  pole  coordinate  of  participant  ZM  (Fig.  4).  Their  values  are 
usually  greater  for  the  x  pole  coordinate  than  for  y  except  participant  VT  (Fig.  3)  and 
ensemble  predictions  (Fig.  5)  where  they  are  of  the  same  order.  A  decrease  of  kurtosis  with 
prediction  length  suggests  that  the  probability  distribution  becomes  more  flat  and  has  larger 
tails  than  a  normal  distribution. 


Days  in  the  future 


Days  in  the  future 


Days  in  the  future 


Days  in  the  future 


- BL-ep  - WK  - ZM  - BL  - VT 


Fig.  6.  Standard  deviations  and  mean  absolute  errors  computed  by  different  participants  of  the 
project  and  of  ensemble  predictions  computed  by  USNO  (BL-ep)  (red). 
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Figure  6  shows  the  mean  prediction  errors  of  pole  coordinate  data  represented  by  the 
standard  deviations  and  mean  absolute  errors  computed  for  different  participants  of  the 
project  and  of  ensemble  predictions  computed  by  USNO.  Note  that  the  mean  prediction  errors 
of  the  ensemble  prediction  are  greater  than  the  smallest  mean  prediction  errors  for  an 
individual  participant  for  the  x  pole  coordinate.  The  mean  ensemble  prediction  error  of  the  y 
pole  coordinate  for  prediction  lengths  less  than  20  days  is  of  the  order  of  the  smallest 
prediction  error  for  an  individual  participant. 

To  explain  the  skewness  variations  for  different  participants,  the  differences  between 
pole  coordinate  data  and  their  least-squares  +  autoregressive  (LS+AR)  predictions  (Kosek  et 
al.  2004)  of  WK  were  computed  (Fig.  7).  In  this  prediction  method  the  IERS  eopc04. 62-now 
data  were  used  (IERS  2011).  The  skewness  for  the  x  pole  coordinate  predictions  is  usually 
less  than  for  y  (Figs.  1-5)  and  these  differences  show  that  the  sign  for  x  is  mostly  negative 
while  the  sign  for  y  is  mostly  positive  during  the  EOPCPPP.  Note  that  between  the  EOPPCC 
and  the  EOPCPPP  the  sign  of  these  differences  for  the  x  and  y  pole  coordinate  is  mostly 
negative  which  indicates  that  the  tail  on  the  left  side  of  the  probability  density  function  is 
longer  than  the  right  side,  and  it  deviates  significantly  from  a  normal  distribution.  Note  that 
between  2004  and  2009  the  change  of  sign  of  these  differences  for  the  y  pole  coordinate 
shows  a  quasi-annual  variation.  Because  the  last  10  years  of  pole  coordinate  data  is  used  in 
the  least-squares  model  of  the  LS+AR  prediction  (Kosek  et  al.  2004),  the  constant  phase  of 
the  annual  oscillation  in  the  extrapolation  model  represents  the  average  phase  for  this  time 
span  of  data.  This  constant  phase  in  the  extrapolation  model  is  not  equal  to  its  real  value  for 
the  prediction  epoch  due  to  significant  phase  variations  of  the  annual  oscillation  in  pole 
coordinate  data  (Kosek  et  al.  2006). 


years 

Fig.  7.  The  differences  between  IERS  pole  coordinate  data  and  their  LS+AR  predictions 

computed  by  WK  for  90  days  in  the  future. 


It  was  found  that  the  short  term  (up  to  100  days  in  the  future)  mean  prediction  errors  of 
pole  coordinate  data  are  of  the  same  order  as  the  mean  prediction  errors  of  the  pole  coordinate 
model  data  computed  from  joint  atmospheric-ocean  excitation  functions  (Kosek  2010,  Kosek 
et  al.  2011).  The  semblance  function  (Kosek  et  al.  2011)  and  polarization  function  (Kosek 
2010)  of  pole  coordinate  data  and  pole  coordinate  model  data  computed  from  the  joint 
atmospheric-ocean  excitation  functions  are  very  similar  especially  for  the  annual  frequency 
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band  which  means  that  the  annual  oscillation  in  pole  coordinate  data  is  mostly  excited  by  the 
joint  atmospheric-ocean  excitation  function.  Thus,  the  main  cause  of  short-term  prediction 
errors  of  pole  coordinate  data  can  be  due  to  irregular  phase  variations  of  the  annual  oscillation 
in  the  atmospheric-ocean  excitation.  The  autoregressive  prediction  of  the  least-squares 
extrapolation  residuals  was  not  able  to  take  these  phase  variations  of  the  annual  oscillation 
into  account. 

5.  CONCLUSIONS 

The  mean  ensemble  prediction  error  of  the  x  pole  coordinate  is  larger  than  the  smallest 
prediction  error  for  an  individual  participant  and  the  mean  ensemble  prediction  error  of  the  y 
pole  coordinate  for  prediction  lengths  less  than  20  days  is  of  the  order  of  the  smallest 
prediction  error  for  individual  participant. 

The  skewness  values  for  different  participants  show  that  the  probability  distribution 
becomes  more  nonsymmetrical  when  the  prediction  length  increases.  Usually,  during  the 
EOPCPPP  the  signs  of  skewness  are  opposite  for  x  and  y  pole  coordinate  predictions. 

The  kurtosis  values  usually  decrease  with  the  prediction  length,  which  means  that  the 
probability  distribution  becomes  more  flat  and  has  larger  tails  than  a  normal  distribution. 

The  increase  of  short-term  prediction  errors  of  pole  coordinate  data  can  be  caused  by 
irregular  variations  of  the  annual  oscillation  phase  in  polar  motion  and  its  atmospheric-ocean 
excitation. 
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